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Abstract 

The well-known Edwards- Wilkinson equation with a flow term added exhibits 
a smoothing fixed point in addition to the normal EW fixed point. Based 
on this, we present a model of sandpiles involving a coupling between fixed 
and mobile grains, which then shows the smoothing behaviour that would be 
expected to obtain after avalanche propagation. 
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The dynamics of sandpiles have intrigued researchers in physics over recent years [jl],|^ 
with a great deal of effort being devoted to the development of techniques involving for 
instance cellular automata 0,^, continuum equations and Monte Carlo schemes ^ 

to investigate this very complex subject. However what have often been lost sight of in all 
this complexity are some of the extremely simple phenomena that are exhibited by granular 
media that still remain unexplained. 

One such phenomenon is that of the smoothing of a sandpile surface after the propagation 
of an avalanche IQ. It is clear what happens physically: an avalanche provides a means of 
shaving off roughness from the surface of a sandpile by transferring grains from bumps to 
available voids [Q, and thus leaves in its wake a smoother surface. However, surprisingly, 
researchers have not to our knowledge come up with a model of sandpiles that has exhibited 
this behaviour. 

In this paper we present therefore a very simple model involving coupled continuum 
equations which exhibits smoothing. In work that is presently under completion we 
have constructed a cellular automaton model of a sandpile surface where this smoothing 
behaviour is also made manifest. 

This paper comprises : (i) an introductory section, where the equations are presented and 
explained (ii) a section where numerical and analytical results are obtained (iii) a concluding 
section where the physics of our results is discussed. 



I. INTRODUCTION: THE EQUATIONS 

Coupled equations have been seen to be immensely helpful in representing the dynamics 
of systems where there are two clear relaxational mechanisms, for instance those correspond- 
ing respectively to particles moving independently of each other and to their collective motion 
within clusters In the case of sandpiles, these mechanisms correspond respectively to 

the intercluster and intracluster motion of individual grains 

The model we present here involves just such a pair of coupled equations, where the 
equation governing the evolution of clusters ("stuck" grains) is closely related to the very 
well-known Edwards- Wilkinson (EW) model |T0|. The equations are: 

dh 

— = DhVh + cVh + r^ix^t) 

^ = DpV'p - cVh (1) 

where the first of the equations describes the height h{x, t) of the sandpile surface at (x, t) 
measured from some mean (h) , and is precisely the EW equation in the presence of the flow 
term cV/i. The second equation describes the evolution of flowing grains, where p{x, t) is 
the local density of such grains at any point {x,t). As usual, the noise ri{x,t) is taken to be 
Gaussian so that: 

{r]{x, t)r]{x', t')) = A^6{x - x')6{t - t'). 

with A the strength of the noise. Here, (■ ■ ■) refers to an average over space as well as over 
noise. 



Before further analysis, we review some general facts about rough interfaces 11 1. Three 



critical exponents, a, /3, and characterise the spatial and temporal scaling behaviour of 
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a rough interface. They are conveniently defined by considering the (connected) two-point 
correlation function of the heights 

5(x - x', t - t') = (/i(x, t)/i(x', t')) - (/i(x, t)){h{^', t')). (2) 

We have 

^(x,0)~|x|2" (|x|^cx)) and S{0,t) r-. \t\^^ ^ oo), (3) 
and more generally 

s{^,t)^\^rF{\t\/\^n (4) 

in the whole long-distance scaling regime (x and t large). The scaling function F is universal 
in the usual sense; a and z = a/[3 are respectively referred to as the roughness exponent 
and the dynamical exponent of the problem. 

In addition, we have for the full structure factor which is the double Fourier transform 

S{k, uj) 

S{k, uj) ~ uj-^k-^-^"<^{uj/k') (5) 
which gives in the limit of small k and u, 

S{k, U = 0)^ f.-l-2a^z ^ g^f^ = Q, w) ~ uJ-^-^f^-^/^ (a; ^ 0) (6) 

The scaling relations for the corresponding single Fourier transforms are 

S{k, t = 0) ~ A;-i-2" {k 0) and S{x = 0, cj) ~ uj~^~^f^ {u 0) (7) 



II. NUMERICAL AND ANALYTICAL RESULTS 

A. Analysis of the decoupled equation: the Edwards- Wilkinson equation with flow 

For the purposes of analysis, we focus on the first of the two coupled equations (Eq.|l]) 
presented above, 

dh 

— = DhV^h + cVh + r]{x,t) (8) 

noting that this equation is essentially decoupled from the second. (This statement is, 
however, not true in reverse, which has implications to be discussed later). We note that 



this is entirely equivalent to the Edwards- Wilkinson equation |T0[ in a frame moving with 
velocity c 

x' = X + ct , t' = t 

and would on these grounds expect to find only the well-known EW exponents a = 0.5 and 
f3 = 0.25 [|l^]. This would be verified by naive single Fourier transform analysis of Eq.^ 
which yields these exponents via Eq.|^. 
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Eq.H can be solved exactly as follows. Using the standard form 

Lh = r]{x, t) 

with 

Z = - D.V^ - cV] 

we find for the propagator G{k,uj) 

G = {-iuj + Dhk'^ + ikc)~^ 

This can be used to evaluate the structure factor 

{h{k,uj)h{k',u')) 
^^^^^^ = 6ik + t)5i. + .') 

which is the Fourier transform of the full correlation function: 

G{x -x',t- t') = {h{x, t)h{x', t')). (10) 
The solution for S{k, u) so obtained is: 

^^^^''^ = i.-ckr+Dik^ ^''^ 

This is illustrated in Fig.|l| while representative graphs for S{k,uj = 0) and S{k = 0,u) 
are presented in Figs.^ and ^ respectively. It is obvious from the above that S{k,uj) does 
not show simple scaling. More explicitly, if we write 



with ko = c/Dh, and ujq = c^/Dh, we see that there are two limiting cases : 

• for A; ^ kf), S^^{k,u = 0) ~ k^; using again S^^{k = 0,ct;) ~ u"^, we obtain a = 1/2 
and /? = 1/4, z = 2 via Eqs.|. 

• for /c -C ko, S~^{k,u! = 0) ~ /c^;using the fact that the limit S^^{k = 0,0;) is always 
tu^, this is consistent with the set of exponents a = 0, /3 = and z = 1 via Eqs.H. 



The first of these contains no surprises, being the normal EW fixed point, [jlO| while the 
second represents a new, 'smoothing' fixed point. 

We now explain this smoothing fixed point with a simple physical picture. The competi- 
tion between the two terms in Eq.|| determines the nature of the fixed point observed: when 
the diffusive term dominates the flow term, the canonical EW fixed point is obtained, in the 
limit of large wavevectors k. On the contrary, when the flow term predominates, the effect 
of diffusion is suppressed by that of a travelling wave whose net result is to penalise large 
slopes; this leads to the smoothing fixed point obtained in the case of small wavevectors k. 
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We add in passing that this crossover would not have been as transparently obvious 
had we been using the single Fourier transforms S{k, t = 0) and S{x = 0, u) for numerical 
purposes. We illustrate this by writing explicitly the expressions for the relevant quantities: 

5(fc,t = 0)~^ (12) 
S(x = 0, u;) ~ — for u small (13) 
S{x = 0, cu) ~ — j-g- for LJ large (14) 

The examination of S{k,t = 0) (Fig.^ on its own yields no indication of the crossover 
to the smoothing fixed point, while only a detailed examination of the S{x = 0,00) graph 
(Fig.^a) reveals this, indicating a crossover from uj~^'^ at high frequencies to a;~^ at low 
frequencies. Clearly, though this is more time-consuming to obtain, the double Fourier 
transform, ie the full structure factor, provides a much more unambiguous picture of the 
crossover in this problem. 

Additionally, the investigation of even the single Fourier transform to detect the crossover 
is not without its complexities. The detection of the smoothing fixed point, characterised 
by 2 = 1 , needs a careful investigation of the low frequency part of S{x = 0,uj). This 
is elaborated on in Appendix ^ where we discuss the anomalous oscillatory behaviour 
obtained to do with the competition between grid sizes and some characteristic lengths in 
this problem (Appendix 

In view of the above, it is preferable to use the double Fourier transform to obtain an 
unambiguous picture of the structure factor although this strategy might on first appearance 
seem to be a computational overkill. The overwhelming advantage is that, by scanning the 
structure factor as a function of frequency 00 for a fixed k, one immediately sets two frequency 
scales ck and D^k"^, thus making it possible to pick up the relevance of these scales in S{k, u). 

Before turning to an interesting physical application of the two fixed points in our linear 
equation, we mention in passing that our discussion is equally applicable to the Kardar- 



Parisi-Zhang (KPZ) equation |jT2| with the addition of a flow term. Here too, the use of the 



double Fourier transform reveals the presence of the 'smoothing' fixed point due to the flow 



term ||13||. 



B. Coupled equations: a model of smoothing 

Returning to the coupled equations (|l]), we realise from the above that the interface h 
is smoothed because of the action of the flow term which penalises the sustenance of finite 
gradients V/i. As we have seen from the above analysis, this is indeed the case for the 
equation in h, which is effectively decoupled from the equation in p. However, the equation 
in p is manifestly coupled to h so we would need to ensure that no instabilities are generated 
in this, in order for the coupled equations to qualify as a valid model of sandpile dynamics. 

In this spirit, we look first at the value of p averaged over the pile, as a function of time 
(Fig.^a). We observe that the incursions of (p) into negative values are limited to relatively 
small values, suggesting that the addition of a constant background of p exceeding this 
negative value would render the coupled system meaningful, at least to a first approximation. 
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In order to ensure that this average over the pile does not involve wild fluctuations over the 
body of the pile, we examine the fluctuations in p, viz.-y/ {p^) — (p)^ (Fig.Hb). The trends 
in that flgure indicate that this quantity appears to saturate, at least upto computationally 
accessible times. Finally we look at the minimum value of p at any point in the pile over 
a large range of times (Fig.||c); this appears to be bounded by a modest (negative) value 
of 'bare' p. Our conclusions are thus that the fluctuations in p saturate at computationally 
accessible times and that the negativity of the fluctuations in p can always be handled by 
starting with a constant po, a constant 'background' of flowing grains, which is more positive 
than the largest negative fluctuation. 

Physically, then, the above implies that at least in the presence of a constant large density 
Po of flowing grains, it is possible to induce the level of smoothing corresponding to the flxed 
point a = /5 = 0. This model is thus one of the simplest possible ways in which one can 
obtain a representation of the smoothing that is frequently observed in experiments on real 
sandpiles after avalanche propagation p|. 

III. CONCLUSIONS 

In conclusion we have presented a very simple model of sandpile dynamics which man- 
ifests the interfacial smoothing commonly observed after avalanche propagation, both in 
experiments [Q and physically in sand dunes. This consists simply of a set of coupled con- 
tinuum equations which model the exchange process between the stuck grains at the 
interface h and the flowing grains p which flow across the interface. The equation in h is 
closely related to the well-known Edwards- Wilkinson equation [|l^], modifying this simply 
by the addition of a flow term which penalises excess gradients across the interface, and 
thereby causes the smoothing. 

We flnd, on analysing the equations, that there are two flxed points which characterise the 
interface equation: one, the normal diffusive flxed point corresponding toa = l/2,/5 = l/4 
and z = 2, and the other, the new smoothing flxed point corresponding to a = 0, /5 = and 
z = 1. We mention here that the flrst of these is the one that would also be observed in the 
frame of the avalanche moving with the flow velocity c, while the smoothing would only be 
observed in the stationary frame of an observer, who would observe the smoother surface 
left in the wake of the avalanche. 

Finally, in work that is in progress, we have observed quantitatively this smoothing 
phenomenon in a cellular automaton model of sandpiles We hope that this will motivate 
actual experiments to measure more quantitatively this so far qualitatively observed and 
intuited behaviour at sandpile surfaces. 
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APPENDIX A: OSCILLATIONS 



We explain here the origins of the observed oscillations in the structure factor S{x = 0,uj) 
The single Fourier transform S{x, u) is defined by 



1 f°° 

S{x = 0,uj) = — / S{k,Lj)exp{ikx)dk 
zn J~oo 

In the limit of small u the integral can be written as: 



S{x = 0, uj) 



1 

2^ 



DP 



6{uj — ck)dk 



(Al) 



(A2) 



This is the origin of the ballistic behaviour of the flow term and is responsible in part for 
some of the observed oscillations, as we explain below. It is clear from Eqs.|l3| and |14| that 
the crossover from small to large u for a; = should not involve any imaginary quantities, 
and therefore strictly speaking we should not see any oscillatory behaviour in the structure 
factor in this limit. However it is important to realise that the full form of the structure 
factor S{x, u) for finite x does contain imaginary portions, in order to understand fully the 
origin of the obtained oscillations. 

The characteristic length and time scales in our problem are given by (Appendix P) 

D , D 

to = — and Xq = — 

c 

Whenever grid sizes in time or space are comparable to these characteristic lengths, the pro- 
file fluctuates across these intervals, which is then aggravated by the shock fronts associated 
with the flow term. This results in: 

• oscillatory behaviour arising from the non-zero intervals in x associated with the sam- 
pling of the profile to generate the Fourier transform, S{x = 0,lj), which introduce a 
flavour of S{x, u) for finite x. 

• oscillations which become increasingly violent as c increases because of the increased 
fluctuations associated with the ballistic flow term over the grids. 

This line of reasoning is borne out by the following table: 

Parameters 



Grid Sizes & Observation 
Characteristic Lengths 



c=l D, = 1 



c=5 Dh = 1 



c=10 Dh =1 



Ax 


= 0.5, Xq = 


= 1.0 


No Oscillation . 


At 


= 0.001, to 


= 1.0 


Fig.||a 


Ax 


= 0.1, xo = 


= 0.2 


Oscillation(l peak). 


At 


= 0.001, to 


= 0.04 


Fig.||b 


Ax 


= 0.1, Xq = 


= 0.1 


Oscillation( 2 peaks) 


At 


= 0.001, to 


= 0.01 


Fig.|c 



7 



Thus in order to avoid these oscillations, one should choose grid sizes Ax and At in such 
a way that they are always less than the characteristic scales in the problem, i.e., 



Ax <^ Xq and At <^ to. 



APPENDIX B: LENGTH AND TIME SCALES OF INTEREST 

The Edwards- Wilkinson equation with a flow term (cV/i) is given by 

dh 

— ^DhV^h + cWh + r){x,t) (Bl) 
ot 

with 

{r]{x, t)r]{x', t')) = A^S{x - x')S{t - t') 
We perform a scale change in order to make this equation dimensionless using 

t — toT , X — xqX , h — hoH 



and obtain 



where. 



^h^^h—, C = C— , do=—\ — 
Xq Xq flo V 



and 



iax, T)^{X', T')) = 5{X - X')5{T - T') 



Next we make a choice of the scales (xq, to and ho) such that the coefficients iP^, c' and 5q 
are unity. This gives the following scales 

= — , xq^ — , ho^ — (B3) 

C C -y C 

An estimate of the exponent z can be made by considering the diffusion and flow term 
separately. Substituting L for xq along with the condition D'j^ = c' = 1 we realise that on 
length scales comparable to the system size, L, the associated time scale for diffusion is 
Li^ /Dh and that for flow is given by L/c. The physical meaning of these time scales is the 
time taken for the diffusive and flow-like correlations to span the system. This is the time 
(tx ~ L^) for crossover to the saturated regime. Clearly the above time scales indicate z = 2 
for diffusion-dominated behaviour and 2; = 1 for flow-dominated behaviour respectively. 
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FIGURES 



FIG. 1. The double Fourier transform, S{k = ki,uj) for the h-h correlation function for three 
wavevectors h = 0.02, k2 = 0.08 and ks = 0.12. (c = 2 and Dh = A'^ = 1.0) 

FIG. 2. The double Fourier transform, S{k,uj = 0) for the h-h correlation function, along with 

a fit to (A:2-00±0.03 + ^4.00±0.05)-l_ ^^^j^^^^2^ 

FIG. 3. The double Fourier transform, S{k = 0, to) for the h-h correlation function, along with 
a fit to u;-2.oo±o.05, (c = D;, = = 1.0) 

FIG. 4. The single Fourier transform S{k, t = 0) for the h-h correlation function, along with 
a power-law fit to /c-i-97±o.02_ 

FIG. 5. (a)The single Fourier transform S{x = 0, to) for the h-h correlation function, along 
with a power-law fit to (_^;-i-9'i'±o.05 £qj. gj^^g^^ ^ g^^id (^-i-48±o.02 large uj. 

(b) and (c) show the oscillatory behaviour (see Appendix ^) of S{x = 0, m) near the crossover 
region. The oscillations increase with increasing c , as can be seen from a comparison of Fig. (b) 
(c = 10) and Fig. (c) (c = 5). Dh = = 1.0 throughout. 

FIG. 6. Curves relating to the behaviour of the equation in p (Eq.^) over 10^ timesteps. 

(a) shows the behaviour of the mean (p) as a function of time, where the average over the sandpile 
was computed over a sample configuration. 

(b) shows the variance Prms of p as a function of time. This was computed for 100 configurations 
over the sandpile surface. 

(c) shows the bounds of the fluctuations in p, pmin and p 

max •) ^ 

function of time. 
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